# Code for Figure H1: Instrument Variation

rm(list = ls())

## ---------------------------------------
## Load Packages 
## ---------------------------------------

library('AER')
library('ivpack')
library('data.table')

## ---------------------------------------
## Load Data and Functions
## ---------------------------------------
load("../data0812.RData")

directory <- "../functions/"
functions <- list.files(directory)  
loadfunctions <- sapply(functions, FUN = function(x)source(paste0(directory, x)))

## ---------------------------------------
# Construct Instrument
## ---------------------------------------
data0812 <- constructIV(data0812)

## ---------------------------------------
#  Select Last Case Before Election
## ---------------------------------------
data0812 <- lastCase(data0812)


## ---------------------------------------
## Two Stage Least Squares Specifications
## ---------------------------------------

time.controls <- "as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow) + as.factor(court_shift) + as.factor(totOGS2)"

outc.1 <- "vote2012"
endo.1 <- "pti"
inst.1 <- "judgeiv"

pdf("../output/Figures/Figure_H1.pdf", h=8,w=16)    
par(mfrow=c(1,2), mar = c(4, 6, 4, 6), mgp = c(2, 1, 0))
form.3 <- formula(paste(endo.1, "~", time.controls))
form.2 <- formula(paste(inst.1, "~", time.controls))

mod1 <- lm(form.3, data = data0812)
data0812$pti_res <- data0812$pti - predict(mod1, data0812)

mod1b <- lm(form.2, data = data0812)
data0812$judgeiv_res <- data0812$judgeiv - predict(mod1b, data0812)

mod2 <- lm(pti_res ~ judgeiv_res + I(judgeiv_res^2), data = data0812)
judgeiv_res <- seq(min(data0812$judgeiv_res), max(data0812$judgeiv_res), by = 0.001)
data.pred2 <- data.table(judgeiv_res)
myPredict <- predict(mod2, newdata = data.pred2, interval = "confidence")

hist(data0812$judgeiv_res, 
     freq = F, breaks = 10, las = 1,
     xlab = "Leave-out-case Residualized Pretrial Detention Rate",
     main = "", ylim = c(0, 16))
x <- judgeiv_res
ix <- sort(x,index.return=T)$ix

par(new = T)
plot(x, myPredict[, 1], axes=F, xlab=NA, ylab=NA, type = "l", col = "red")
polygon(c(rev(x[ix]), x[ix]), c(rev(myPredict[ ix,3]), myPredict[ ix,2]), col = rgb(0.7,0.7,0.7,0.4) , border = NA)
lines(x[ix], myPredict[ix , 1], col=2, lwd=2)
axis(4, at = seq(from = min(myPredict), to = max(myPredict), by = 0.01), 
     labels = round(seq(from = min(myPredict), to = max(myPredict), by = 0.01), 2), 
     col = 'red', col.axis = 'red', las = 1)
mtext("Residualized Pretrial Detention Rate", side = 4, col = "red", line = 4)

form.3 <- formula(paste(outc.1, "~", time.controls))
mod1 <- lm(form.3, data = data0812)
data0812$t_res <- data0812$vote2012 - predict(mod1, data0812)

mod2 <- lm(t_res ~ judgeiv_res + I(judgeiv_res^2), data = data0812)

judgeiv_res <- seq(min(data0812$judgeiv_res), max(data0812$judgeiv_res), by = 0.001)
data.pred2 <- data.table(judgeiv_res)
myPredict <- predict(mod2, newdata = data.pred2, interval = "confidence")

hist(data0812$judgeiv_res, 
     freq = F, breaks = 10, las = 1,
     xlab = "Leave-out-case Residualized Pretrial Detention Rate",
     main = "", ylim = c(0, 16))
x <- judgeiv_res
ix <- sort(x,index.return=T)$ix

par(new = T)
plot(x, myPredict[, 1], axes=F, xlab = NA, ylab = NA, type = "l", col = "red", ylim = c(-0.02, 0.02))
polygon(c(rev(x[ix]), x[ix]), c(rev(myPredict[ ix,3]), myPredict[ ix,2]), col = rgb(0.7,0.7,0.7,0.4) , border = NA)
lines(x[ix], myPredict[ix , 1], col=2, lwd=2)
axis(4, at = seq(from = min(myPredict), to = (max(myPredict)), by = 0.01), 
     labels = round(seq(from = min(myPredict), to = (max(myPredict)), by = 0.01), 2), 
     col = 'red', col.axis = 'red', las = 1)
mtext("Residualized Turnout", side = 4, col = "red", line = 4)

dev.off()

cat("\nFigure H1: Instrument Variation\n")
cat("Saved Figure H1 in /output/Figures/Figure_H1.pdf\n\n")
